Sprach und Signalverarbeitung

Für diese Aufgabe verwenden wir die Pakete WAV um Klangdateien auszulesen und Plots für das Plotten. Zudem verwenden wir plotly für interaktive Plots.

In [54]:
#Pkg.add("WAV")
#Pkg.add("Plots")
using WAV,Plots

plotly()
Out[54]:
Plots.PlotlyBackend()

Im Folgenden arbeiten wir mit zwei Klangdateien.

In [55]:
(s1,fs1) = wavread("../speech1.wav");
(s2,fs2) = wavread("../speech2.wav");

Aufgabe 1

(a)

Die Abtastrate ist der zweite Ausgabeparameter der Funktion wavread

In [56]:
println("Abtastrate für speech1.wav: ",fs1," Hz")
println("Abtastrate für speech2.wav: ",fs2," Hz")

Abtastrate für speech1.wav: 16000.0 Hz
Abtastrate für speech2.wav: 16000.0 Hz

(b)

If s is a signal with sampling frequency fs, then the i-th measurement in s takes place at time (i-1)*period_length where period_length = 1/fs.

(Note: Indexing of vectors starts at 1 in both Julia and Matlab)

In [57]:
v_sample_time1 = (0:(length(s1)-1))/fs1;
v_sample_time2 = (0:(length(s2)-1))/fs2;
In [58]:
plot(v_sample_time1,s1;title="speech1.wav",xlab="Zeit [s]")
Out[58]:
In [59]:
plot(v_sample_time2,s2;title="speech2.wav",xlab="Zeit [s]")
Out[59]:
  • Die stillen Regionen zeichnen sich durch eine Amplitude nahe 0 aus.
  • stimmhaft: Sehr periodisch, relativ hohe Amplitude
  • stimmlos: Hohe Unregelmäßigkeit (bei jeder Skalierung), Amplitudenmaxima variieren mehr

TODO

(c)

  • Periode von 0.3428 bis 0.3475 (Schätzung)
  • Frequenz: 1/(0.3475 - 0.3428) ~ 212 Hz

TODO

Aufgabe 2

Gegeben sind

  • ein Signal v_signal
  • seine Abtastrate sampling_rate
  • die Länge eines Zeitfensters frame_length
  • Abstand zwischen den Startpunkten zweier aufeinanderfolgender Zeitfenster frame_shift (uff, das geht sicher eloquenter)

Die erste Aufgabe besteht darin die Abmessungen für unsere Zeitfenster, die in Sekunden angegeben sind, als eine Anzahl von Samples auszudrücken. Dadurch ergeben sich folgende Definitionen:

  • samples_per_frame = Int(floor(frame_length * sampling_rate))
  • samples_per_shift = Int(floor(frame_shift * sampling_rate))

Als nächstes stellt sich die Frage wie viele Fenster in unser Signal passen. Wir nehmen sinnvollerweise an, dass in unser Signal mindestens ein Fenster (genauer gesagt, das erste Fenster) passt. Jedes darauf folgende Fenster konsumiert von den restlichen length(v_signal) - samples_per_frame Samples genau samples_per_shift. Dementsprechend ergibt sich die Formel

  • shift_count = 1 + Int(floor((length(v_signal) - samples_per_frame)/samples_per_shift))

Diese Daten reichen nun aus, um die Ausgabematrix zu befüllen.

Für die Mittelpunkte der Zeitfenster ergibt sich durch folgende Überlegung: Der Mittelpunkt des ersten Fensters befindet sich bei frame_length/2. Jeder darauffolgende Mittelpunkt ergibt sich aus einer Verschiebung um frame_shift.

  • v_time_frame = frame_length/2 + frame_shift * (0:(shift_count-1))
In [60]:
my_windowing = function(v_signal, sampling_rate, frame_length, frame_shift)
    samples_per_frame = Int(floor(frame_length * sampling_rate))
    samples_per_shift = Int(floor(frame_shift * sampling_rate))
    shift_count = Int(floor((length(v_signal) - (samples_per_frame - samples_per_shift))/samples_per_shift))

    m_frames = zeros(samples_per_shift, shift_count)
    for j in 1:shift_count
        for i in 1:samples_per_shift
            m_frames[i,j] = v_signal[samples_per_shift*(j-1)+i]
        end
    end
    v_time_frame = frame_length/2 + frame_shift * (0:(shift_count-1))

    (m_frames, v_time_frame)
end
Out[60]:
(::#29) (generic function with 1 method)

Aufgabe 3

(a)

In [61]:
(frames1,frame_times1) = my_windowing(s1,fs1,32e-3,16e-3);
(frames2,frame_times2) = my_windowing(s2,fs2,32e-3,16e-3);

(b)

In [62]:
autocor = x -> conv(x,reverse(x))/length(x)
Out[62]:
(::#31) (generic function with 1 method)

(c)

truncated version of autocor

In [63]:
autocor_t = x -> autocor(x)[1:length(x)]
Out[63]:
(::#33) (generic function with 1 method)

(d)

In [64]:
ff_estimator = function(s,fs)
    (frames, frame_times) = my_windowing(s,fs,32e-3,16e-3)
    (samples_per_frame, shift_count) = size(frames)
    
    ac = zeros(frames)      # auto correlations
    freq = zeros(frame_times) # ff estimate per frame
    
    # 80Hz:400Hz frequency window
    lb = Int(floor(fs/400))  # lower bound
    ub = Int(floor(fs/80))   # upper bound
    
    for j in 1:shift_count
        ac[:,j] = autocor_t(frames[:,j])
        period = lb - 1 + indmax(ac[lb:ub,j])   # period estimate (measured in number of samples)
        freq[j] = fs/period
    end
    
    return (frame_times, freq)
end
Out[64]:
(::#35) (generic function with 1 method)

(e)

Signal strengths s1 and fundamental frequency estimates ff1 have different units (and magnitudes) which is why we have to rescale in order to see the details of both data in one plot.

In [68]:
(tf1,ff1) = ff_estimator(s1,fs1);
plot([v_sample_time1,tf1],[s1*2000,ff1])
Out[68]:
In [67]:
(tf2,ff2) = ff_estimator(s2,fs2);
plot([v_sample_time2,tf2],[s2*2000,ff2])
Out[67]:
  • The estimator gives reasonable frequencies for voiced regions and much less so in unvoiced regions, where it fluctuates a lot.
  • Our manually determined estimate and the computed estimate differ by a factor of 2 (!!!!) :(

TODO